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Abstract 



Understanding the differentiation, a biological process from a multipo- 
tent stem or progenitor state to a mature cell is critically important. We 
develop a theoretical framework to quantify the underlying potential land- 
scape and biological paths for cell development and differentiation. We 
propose a new mechanism of differentiation and development through bind- 
ing/unbinding of regulatory proteins to the gene promoters. We found in- 
deed the differentiated states can emerge from the slow promoter binding/ 
unbinding processes. Furthermore, under slow promoter binding/unbinding, 
we found multiple meta-stable differentiated states. This can explain the ori- 
gin of multiple states observed in the recent experiments. In addition, the 
kinetic time quantified by mean first passage transition time for the dif- 
ferentiation and reprogramming strongly depends on the time scale of the 
promoter binding/unbinding processes. We discovered an optimal speed for 
differentiation for certain binding/unbinding rates of regulatory proteins to 
promoters. More experiments in the future might be able to tell if cells differ- 
entiate at at that optimal speed. In addition, we quantify kinetic pathways 
for the differentiation and reprogramming. We found that they are irre- 
versible. This captures the non-equilibrium dynamics in multipotent stem 
or progenitor cells. Such inherent time-asymmetry as a result of irreversibil- 
ity of state transition pathways as shown may provide the origin of time 
arrow for cell development. 



1 Introduction 



During cell differentiation, the cell evolves from undifferentiated phenotypes 
in a multipotent stem or progenitor state to differentiated phenotypes in a 
mature cell. In this process, the gene regulatory network, which governs the 
progressive changes of gene expression patterns of the cell, forces the cell to 
adopt the cell type-specific phenotypes. Cells can have states with the higher 
probability of appearance, which leads to different cell phenotypes. Different 
cell phenotypes correspond to different basins of attractions on the potential 
landscape (1~3). Therefore the differentiation and developmental process 
of the cell can be thought as the evolution of the underlying landscape 
topography from one basin to to another. One grand challenge is to explain 
how this occurs, what is the underlying mechanism and how to quantify the 
differentiation and developmental process. Furthermore, the unidirectional 
developmental process posses another challenge to explain the origin of time 
arrow. 

In the cell, intrinsic fluctuations are unavoidable due to the limited num- 
ber of protein molecules. There have been increasing numbers of studies on 
how the gene regulatory networks can be stable and functional under such 
highly fluctuating environments (4-6). In addition, the gene state fluctua- 
tions from the regulatory proteins binding/unbinding to the promoters can 
be significant for gene expression dynamics. Conventionally, it was often 
assumed that the binding/unbinding is significantly faster than the synthe- 
sis and degradation (adiabatic limit) (7, 8). This assumption may hold in 
some prokaryotic cells in certain conditions, in general there is no guaran- 
tee it is true. In fact, one expects in eukaryotic cells and some prokaryotic 
cells, binding/unbinding can be comparable or even slower than the cor- 
responding synthesis and degradation (non-adiabatic limit). This can lead 
to nontrivial stable states and coherent oscillations appearing as a result 
of new time scales introduced due to the non-adiabaticity (9-18). There- 
fore, the challenge for us is to understand how the biological differentiation 
and reprogramming can be functional under both intrinsic fluctuations and 
non-adiabatic fluctuations. 

Previous studies showed that the change in the self activation regulatory 
strengths can cause the differentiation of phenotypes (2, 3, 19). In this ar- 
ticle, we used a canonical gene regulatory circuit module to study cell fate 
decision and commitment in multipotent stem or progenitor cells (2, 3, 19). 
We will study a model of cell developmental circuit (Fig. 1) (20) which is 
composed of a pair of mutually inhibiting but self activating genes. This gene 
regulatory motif has been found in various tissues where a pluri/multipotent 
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stem cell has to undergo a binary cell fate decision (21, 22). For example, 
in the multipotent common myeloic progenitor cell (CMP) facing the bi- 
nary cell fate decision between the myeloid and the erythroid fate, the fate 
determining transcription factors (TF), PU.l, and GATAl, which promote 
the myeloid or the erythroid fates, respectively, form such a gene network 
circuit. The relative expression levels A (PU.l) and B (GATAl) of these 
two reciprocal TFs can bias the decision toward either lineage (20, 22). 

We found that the change in the time scale of the binding/unbinding of 
regulatory proteins to the promoters may provide an new important mech- 
anism for the cell differentiation. We studied the underlying potential land- 
scapes associated with the differentiation and developmental process and 
found that the underlying landscapes developed from un-differentiated mul- 
tipotent state to the differentiated states as the binding/unbinding rate de- 
creased to the slow non-adiabatic binding region. In addition, in the slow 
non-adiabatic binding region, we predicted the emergence of multiple meta- 
stable states in the development of multipotent stem cells and explained the 
origin of this observation in the experiments (19). We also calculated the 
mean first passage transition time for the differentiation and reprogramming. 
We found that the mean first passage transition time strongly depends on 
the time scale of the promoter binding/unbinding processes. There is an 
optimal speed for differentiation and development with certain promoter 
binding/unbinding rates. It will be natural to ask whether the differentia- 
tion and development happens at this optimal speed? Future experimental 
and bioinformatics studies might be able to give the answer. We quantified 
the kinetic pathways for the differentiation and reprogramming. We found 
that they are irreversible. This captures the non-equilibrium prosperities for 
the biological processes of the underlying gene regulatory networks in mul- 
tipotent stem or progenitor cells. It may provide the origin of time arrow 
for development. 

2 Methods and Materials 

As shown in Fig. 1, the gene regulatory circuit that governs binary cell fate 
decision module consists of mutual regulation of two opposing fate determin- 
ing master TF A and B. The module has been shown to control develop- 
mental cell fate decision and commitment in several instances of multipotent 
stem or progenitor cells that faces a binary fate decision, (i.e., GATAl and 
PU.l) (21, 22). A and B are coexpressed in the multipotent undecided cell 
and committed to either one of the two alternative lineages is associated 
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with either one factor dominating over the others, leads to expression pat- 
terns in a mutually exclusive manner (21, 23). Importantly, in many cases 
the genes A and B also self-activate (positive autoregulate) themselves (Fig. 
1). Here, the hybrid promoter a can be bound by the regulatory protein 
/3 with the binding rate hai3 and dissociation rate fap (both ha/s and /q/? 
can depend on protein concentration n^). The synthesis of protein a is con- 
trolled by the gene state of promoter a. There are two types of genes, A 
and B, to be translated into proteins A and B respectively. The proteins 
A{B) can bind to the promoter of the gene A{B) to activate the synthesis 
rate of A{B), which makes a self-activation feedback loop. The proteins 
A{B) can bind to the gene B{A) to repress the synthesis rate of B{A), 
which makes a mutual repression loop. Here, both protein A and protein B 
bind on promoters as a dimer with the binding rate ^haAnAin-A — 1) and 
^haB'iT'BiiT'B — 1) respectively. Therefore, each gene has 4 states with self 
activator binding or non-binding and with mutual repression from another 
gene binding or non-binding (assuming we have two different binding sites, 
one for self activator and one for the other gene). The whole system has 
16 gene states in total. For simplicity, we neglect the roles of mRNAs by 
assuming translation processes are very fast. The model can be expressed 
by the following chemical reactions: 

Ol^ + 2A^Ol\ 010 + 2^^000 (1) 

faA faA 

Oi' + 2B^Oi', OOi^2i3^0OO (2) 

faB faB 

O'i^A, 0%^B, A^% B^% (3) 

with a = A{B) for the hybrid promoter of gene A[B). For the gene state 
index ij of gene O^, the first index i = 1(0) stands for the activator protein 
A unbound(bound) on the promoter a; the second index j = 1(0) stands 
for the repressor protein R unbound(bound) on the promoter a. (7^ {9b) is 
the synthesis rate of the protein A (B) when the gene A (B) is in state ij. 
The probability distribution of the microstate is indicated as PijkiinA^nB) 
where ua and ub are the concentration of the activator A and the repressor 
B respectively. The index represents the gene A occupation state by the 
protein A{B) and the index k{l) represents the gene B occupation state by 
the protein A[B). This results sixteen master equations for the probability 
distribution which are shown explicitly in Supporting Material (SM). 

The steady state probability distribution satisfies — = for all 
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k, I. The total probability distribution is = ^^j^i Pj^j^i ■ The gener- 
alized potential function U of the non-equilibrium network can be quantified 
as: U{nA,nB) = — InP^*'*^ It maps to the potential landscape, which gives 
a quantitative measure of the global stability and function of the underlying 
network (25). Above equations are difficult to deal with because each one 
actually represents an infinite number of equations (n range from to oo). 
A direct way to find the steady state P^^ is through kinetic simulations (24). 
Here, we will use Monte Carlo Simulation to find the stead state distribution 
of master equations (see Supporting Material (SM)). 

3 Results 

In our calculations, we only consider the A-B symmetric case: 

hAA = hsB = hA, hAB = hsA = hn (4) 
Jaa = fsB = Ja, Iab = fsA = Ir (5) 
kA = ks = k (6) 

We define the normalized binding/unbinding rate of the gene states: uja = 
fA/k, OJR = fn/k, and equilibrium constants: = fA/fiA, = fn/hR, 
which indicate the ratio between unbinding and binding speed. There are 
four gene states for each gene and the synthesize rates from gene A and B 
are also symmetric: gfj = g^^. When gene A is bound by protein A (self 
activation) while not bound by protein B (mutual repression) , the synthesize 
rate of protein for protein A is the largest: g^^ = Fa + gfi = Fr + g^^ = 
Fa + Fr + g^Q , where Fa is the activation strength and Fr is the repression 
strength. Here, we choose equilibrium constants = X^ = 45, symmetric 
binding/unbinding speed uja = = w, the repression strength Fr = 60 
and scale the time to make k = 1. 

3.1 The potential Landscapes and two mechanisms for cell 
fate decision of development and differentiation 

Such circuits with above control parameters can generate asymmetric attrac- 
tors representing the differentiated states with almost mutually excluding 
expression of protein A (i.e. GATAl) and B (i.e. PU.l). In addition, cen- 
tral symmetric attractor states characterized by approximately equal levels 
of UA and hb expression can also be generated, which represent the mul- 
tipotent state that exhibits the characteristic balanced or promiscuous ex- 
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pression of the two opposing, fate-determining concentrations-a hallmark of 
the indeterminacy of the undecided multipotent stem cell. 

We plotted the potential landscape in ua-ub plane for different activa- 
tion strength Fa and binding/unbinding speed u in Fig. 2(a), 2(b), 2(c), 
2(d), 2(e), 2(f), 2(g), 2(h), 2(i) for contour view, and Fig. 3(a), 3(b), 3(c), 
3(d), 3(e), 3(f), 3(g), 3(h), 3(i) for 3 dimensional view. In these figures, we 
found two kinds of mechanisms for the cell differentiation. 

During the developmental process, the self activation regulation coming 
from an effective regulation and its change is due to the regulations on these 
transcription factors mediated by other regulators such as Klf4. When the 
self activation is strong (large Fa), the system is mono-stable with one 
un-differentiated central basin, as in Fig. 2(a) (or 3(a)). As self activation 
strength Fa decreases, the central basin gets weaker and differentiated basins 
on both sides start to develop, which results tri-stability as in Fig. 2(d) (or 
3(d)). When self activation strength Fa — > 0, the circuit will reduce to 
a normal symmetric toggle switch. For toggle switch, ua and ub can not 
be both large in adiabatic limit, because they suppress each other. Then, 
the un-differentiated central basin disappeared and two differentiated basins 
on both sides survives, which gives bi-stability as in Fig. 2(g) (or 3(g)). 
Therefore, decreasing the self activation regulatory strength Fa will lead 
the cell system to differentiate. Changing of the effective self activation 
regulatory strengths of transcription factors binding to the genes therefore 
provides a possible differentiation mechanism which is currently under study 
(2, 3, 20, 22). 

We would like to point out that there is another possible mechanism 
of the cell differentiation from the slow binding/unbinding of protein reg- 
ulators to gene promoters. We noticed that for a fixed activation strength 
Fa, cells can develop more stable differentiated states on both sides. As 
shown in Fig. 2(a) (or 3(a)), 2(b) (or 3(b)), 2(c) (or 3(c)) and Fig. 2(d) 
(or 3(d)), 2(c) (or 3(e)), 2(f) (or 3(f)), when binding/unbinding rate uj de- 
creases, the un-differentiated central basin becomes weaker and less stable, 
while differentiated basins on both sides become stronger and more stable. 
We also noticed that in the non-adiabatic slow binding limit (small bind- 
ing/unbinding rate co), multiple meta-stable basins show up. In addition, 
in the non-adiabatic slow binding limit, cells have chances to extinct and 
there are "extinct basins" near [ha = 0,nB = 0), as shown in Fig. 2(c) 
(or 3(c)), 2(f) (or 3(f)), 2(i) (or 3(i)). These behaviors are directly due 
to the non-adiabatic effect: slow binding/unbinding of protein regulators 
to promoters. When the binding/unbinding rate u is small, the interac- 
tions (either repressions or activations) between gene states are weak and 
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different gene states statistically co-exist in cells. Each gene state will give 
a basin in the concentration and sum of these basins will lead a multiple 
stable potential landscape. This results to the development and differentia- 
tion with slow binding from the original undifferentiated equally populated 
single basin of attraction with fast binding. Slow binding provides another 
possible mechanism for differentiation and development. 

3.2 Kinetic and optimal speed for development and differ- 
entiation 

To quantitatively characterize the dynamics of the differentiation and the 
reverse process as reprogramming, we study the speed of differentiation and 
reprogramming in terms of mean first passage time (MFPT) , as shown in 
Fig. 4(a), 4(b) and 4(c). In an attractor landscape, the lifetime of an at- 
tractor reflects its stability, which can be measured by MFPT. MFPT is 
the average transition time induced by intrinsic statistical fluctuations of 
molecule numbers between attractors on a landscape, since the traversing 
time represents how easy to switch from one place to another. When the 
binding/unbinding rate u is relatively large, the un-differentiated central 
basin becomes more stable, as in Fig. 2(a) (or 3(a)), 2(d) (or 3(d)), 2(g) 
(or 3(g)), and cells have more chances to stay in the un-differentiated state. 
Therefore, the differentiation process will be more difficult and MFPT is 
longer for faster binding. For the differentiation process, it is noticed that, 
as the binding/unbinding rate oj increases, the MFPT decreases first, and 
then increases. In the non-adiabatic limit (small binding/unbinding rate uj), 
the differentiation limiting step is the binding/unbinding events. Therefore, 
increasing binding/unbinding speed uj will accelerate the kinetics from the 
un-differentiated central basin to differentiated side basins. So for the differ- 
entiation process, caused from faster binding to slower binding of regulatory 
proteins to the genes, we notice that the speed for differentiation is slow 
when state is dominated by undifferentiated state for faster binding, and is 
also slower for slower binding which is due to the occasional binding being 
the rate limiting step for differentiation. There is an optimal speed for dif- 
ferentiation. As binding becomes faster from low speed end (non-adiabatic 
limit), the speed of differentiation is controlled by the binding speed and 
therefore increases. As the binding becomes even faster, the differentiation 
is dominated by the escape from the undifferentiated basin of attraction 
and therefore is significantly slowed down. This creates an optimal speed 
for differentiation and development. 

The reverse process of cell differentiation is the reprogramming of differ- 
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entiated cells back to a multi- or pluripotent state. In Fig. 4(a), 4(b) and 
4(c), the MFPT for the reprogramming for different self activation strength 
Fa and binding/unbinding speed w is plotted. We observed that, for a typi- 
cal differentiated system, as in Fig. 2(c) (or 3(c)) and Fig. 2(g) (or 3(g)), the 
reprogramming chance is very low and requires a very long time. For self ac- 
tivation strength Fa = 20 (Fig. 4(a)) and Fa = 13 (4(b)), the MFPT for the 
reprogramming decrease as the increasing of thbinding/unbinding speed u, 
because the stability of un-differentiated symmetric central state increases 
with the binding/unbinding speed w as we can see in potential landscapes, 
Fig. 2(a) (or 3(a)), 2(b) (or 3(b)), 2(c) (or 3(c)), 2(d) (or 3(d)), 2(e) (or 
3(e)) and 2(f) (or 3(f)). While in Fig. 4(c), since there is no self-activation 
and no stable symmetric central basin in the landscape, the reprogramming 
is difficult and the MFPT is very long for different the binding/unbinding 
speed uj. 

3.3 Biological dynamic pathways of differentiation and re- 
programming 

Both the differentiation and reprogramming can be caused by the change 
of gene regulations during the developmental process. Here we consider the 
evolution of the binding/unbinding rate from fast to slow uj{t) = lOOOe"''* -|- 
0.001 from 1000 — > 0.001 for the differentiation and the evolution of the bind- 
ing/unbinding rate uj{t) = 1000[1 - e"''*] + 0.001 from 0.001 1000 for the 
reprogramming from slow to fast. The transition paths from Gillespie simu- 
lation are plotted in Fig. 5, accompanied with the potential landscapes for 
the binding/unbinding speed uj = 0.001, 1, 1000. It is interesting to observe 
that the biological dynamic paths are irreversible, i.e. the differentiation 
path and reprogramming path are totally different. In the differentiation 
process, the system stay on the multipotent undifferentiated state for a 
while until binding becomes slower. As the binding becomes slower, the un- 
differentiated state becomes less stable. Furthermore, the gene state can be 
switched through binding/unbinding event of regulatory proteins to the pro- 
moters and the system will then be evolved from the undifferentiated basin 
to the differentiated basin of attraction. In the reprogramming process, the 
system will be gradually attracted into the undifferentiated basin as the in- 
creasing of the binding/unbinding rate to. The paths of differentiation do 
not follow the gradient steepest descent of the potential landscape. They 
do not follow the paths of the reprogramming (the reverse differentiation 
process). This irreversibility reflects the underlying non-equilibrium nature 
of the differentiation and developmental network systems (3). It can give us 
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the fundamental understanding of the biological origin of time arrow in cell 
development. 

4 Conclusion 

We developed a theoretical framework to quantify the potential landscape 
and biological paths for cell development and differentiation. We found a 
new mechanism for differentiation. The differentiated state can emerge from 
the slow promoter binding/unbinding processes. We found under slow pro- 
moter binding, there can be many meta-stable differentiated states. This has 
been observed experimentally (19). Our theory gives a possible explanation 
for the origins of those meta-stable states in the experiments. 

We show that the developmental process can be quantitatively described 
and uncovered by the biological paths on the potential landscape and the 
dynamics of the developmental process is controlled by a combination of the 
intrinsic fluctuations of protein concentrations and gene state fluctuations 
through promoter binding. We also show that the biological paths of the re- 
verse differentiation process or reprogramming are irreversible and different 
from the ones of the differentiation process. 

We explored the kinetic speed for differentiation. We found that the 
cell differentiation and reprogramming dynamics strongly depends on the 
binding/unbinding rate of the regulatory proteins to the gene promoters. 
We found an optimal speed for differentiation and development with cer- 
tain binding/unbinding rates of regulatory proteins to the gene promoters. 
An interesting question we may ask is that is the differentiation and devel- 
opment at optimal speed? More experimental and bioinformatics studies 
might be able to pin down the answer. Furthermore, the irreversibility in 
cell development gives biological examples, which can be easily observed in 
experiments, for the understanding of the origin of time arrow in general 
non-equilibrium systems. 
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Supplementary Informations: Master Equations 



16 Master equations for the canonical gene regulatory circuit of two 
mutually opposing proteins that positively self-regulate themselves, as in 
Fig. 1, are given as following: 

dPiui{nA,nB) _ 
dt ~ 

hAA 



[nA{nA - l)]Piiii{nA,nB) + fAAPmii{nA - 2,nB) 
[ns^nB - l)]PiiiiinA,nB) + fABPimi{nA,nB - 2) 
[nA{nA - l)\Piiii{nA,nB) + fBAPiioi{nA - 2,nB) 
[uBiriB - l)]PiiiiinA,nB) + fBBPiiioinA,nB - 2) 



2 

hAB 

2 

hBA 
2 

hsB 
2 

+kA[{nA + l)Pnn{nA + 1,^^) - nA-Piiii(n^, n^)] 
+kB[{nB + l)Puii{nA,nB + 1) - nBPiui{nA,nB)] 
+gu[Piiii{nA - I^ub) - PiiiiinA,nB)] 
+gii[Piiii{nA,nB - 1) - -Pllll(?^A,?^s)] (SI) 



dPioii{nA,nB) _ 
dt ~ 

--^[nyi(nA - l)]Ploll(r^A,?^B) + /aa-Pooh - 2,^^) 
+ 2)(nB + l)]Puu{nA,nB + 2) - /ab-Ploii("-a, ^^b) 
[uAiriA - l)]Pioiiin A, nB) + /ba-Piooi("-/1 - 2, jib) 



2 



[uBinB - l)]Pion{nA,nB) + /bb-Pioio(?M, "-B - 2) 
+A;A[(nA + l)Pioii(nA + 1,"-^) - 71^^1011(72^,715)] 
+A:b[(?ib + l)Piou{nA,nB + 1) - 715^1011(71^, tib)] 
+510 [-P1011 (77-A - 1,^b) - -Ploii(7ia,71b)] 
+fl'fi[Pioii (nA,nB - 1) - Pioii{nA,nB)] (S2) 
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dPmii{nA,nB) 
dt 



^^[(n^ + 2){nA + l)]Piiii (nA + 2,?!^) - /aa-PoiiiI^a, 



+ 2 



^^'^ - l)]Poill(nA,nB) + /aB-P001i("'A,?^B - 2) 



2 

^BA 

2 



[nA(nA - l)]PQiii{nA,nB) + /ba-Poioi("'A - 2,?ib) 

— - l)\Pmii{nA,nB) + fBBPono{nA,nB - 2) 
+/cA[(nA + l)Poiii("'A + l,'"'^) - nAPoiii{nA,nB)] 
+kB[{nB + l)Poiii{nA,nB + 1) - raB-Poiii("-A, "-b)] 
+5oi[-Poiii(f^A - - Poiii{nA,nB)] 

+gu[Poui{nA,nB - 1) - -Poiii("-a,"-b)] 



dPooii{nA,nB) _ 
dt ~ 

+ 2)(nA + l)]Pioii(nA + 2,725) - /aa-Pooii("-a,?t-b) 
+-^[{i^B + '2){nB + l)]Polll("'A,?^B + 2) - /ab-Pooii("-a,"-b) 
[uAinA - 1)]Pooii("'A,"-b) + /ba-Poooi(?^a - 2,nB) 



2 



[uBinB - l)]Pooii{nA,nB) + /bb-Pooio(?^a, "-B - 2) 



2 

+^'a[(?^A + 1)^0011 (f^A + l,f^B) - ?^A-fooll(^^A,^^B)] 
+kB[{nB + l)Pooii (f^A,f^B + 1) - ?^B-Pooll(^^A,^^B)] 
+5'oo[-fooii(?^A - 1,"-^) - PooiiinA,nB)] 
+gii[PooiiinA,nB - 1) - -Pooii(?^a,?^b)] 
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+ 



dPiuo{nA,nB) _ 



dt 



--^[nA{nA - l)]Piiw{nA,nB) + /AA-Pollo(^^A - 2,725) 
■^^[725(725 - l)]Plllo(?^A,?^B) + fABPww{nA,nB - 2) 
■-^[nA(7iA- l)]Plllo(?^A,?^B) + /ba-Piioo("'A -2,71^) 

— [(725 + 2)(7is + l)]Piiii(7iA,7iB + 2) - fBBPnio{nA,nB) 

+kA[inA + l)Piiio(nA + - nAPiiioinA,nB)] 

+kB[{nB + l)Piiio(nA,nB + 1) - nBPiiio{nA,nB)] 
+gii[Piuo{nA - l,nB) - Piiio{nA,nB)] 
+9io[Piiioi'^A,nB - 1) - Piiio{nA,nB)] 



dPioio{nA,nB) _ 
dt ~ 

--^[^A{nA - 1)]Pioio("-a,?t-b) + /aa-PooioI^-a - 2,71b) 
+ 2)(?2B + l)]Plllo(7^A,?^B + 2) - /ab-Pioio(?M, ?^b) 

--y^["-A("'A - 1)]Pi010("'A,"-b) + /bA-Pi000("'A - 2,71b) 

H — ^[("^ + 2)(7iB + l)]Ploll(7^A,?^B + 2) - /bb-Pioio(?M, ?^b) 
+^a[("'A + l)-Pioio(f^A + l,f^B) - nAPioioinAjns)] 
+kB[inB + l)Pioio(f^A,'^B + 1) - nBPioio{nA,nB)] 
+9io[Pwwin-A -I^ub) - PioioinA,nB)] 
+9io[Pwioi'^A,nB - 1) - Pioio{nA,nB)] 
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dPouo{nA,nB) _ 
dt ~ 

+ + 2){nA + l)]Piiio(raA + 2,nB) - /Ayl-Poiio("-yl, "-b) 

--^[n-BinB - l)]Pollo(^^A,?^B) + /AB-Poolo(?^A, ""-b - 2) 

--^[nAin-A - l)]Poiioin A, ns) + fBAPoiooinA -^^ub) 

— ^[("--s + 2)(?^B + l)]Poin{nA,nB + 2) - /bbPoiio("-a, 

+fe^[(nA + l)Poiio{nA + l,nB) - nAPoiioinA,nB)] 
+kB[{nB + l)Poiio(f^A,'^B + 1) - nBPoiio{nA,nB)] 
+goi[Pouo{nA - l,nB) - Poiio{nA,nB)] 
+9io[Poiw{n-A,nB - 1) - Poiio{nA,nB)] (S7) 



dPooio{nA,nB) _ 
dt ~ 

+-^[{iT'A + 2)(n^ + l)]Pioio(?T-A + 2,nB) - /aa-Pooio("'A, "-s) 
+ 2)(nB + l)]Pollo("'A,?^B + 2) - /ab-Pooio(?M, 

- 1)]Pooio("'A,"-b) + fsAPooooinA - 2,nB) 

^BB 

H — ^[("^ + 2)(nB + l)]Poon{nA,nB + 2) - /bbPooio(?M, 
+/i;^[(n^ + l)Pooio{nA + - f^A-fooio(?^A, 

+kB[inB + l)Pooio(f^A,ns + 1) - nBPooio{nA,nB)] 
+5'oo[-Pooio("'A - - -Poolo(^^A,?^B)] 

+5'fo[-Pooio("-A,raB - 1) - Pooio{nA,nB)] (S8) 
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dPiioi{nA,nB) 



dt 

hAA 



2 

hAB 
2 



[uAiriA - l)]Puoi{nA,nB) + /aa^'oioi - 2,71^) 
[nBiriB - l)]PiioiinA,nB) + /abAooi(?M, - 2) 



+ 2)(nA + l)]Piiii(nA + 2,nB) - /ba-Piioi("'A, 

-^[nB{nB - l)]Piim{nA,nB) + /bb-Piioo(?M, - 2) 

+ l)Piim{nA + ^,nB) - nAPiim{nA,nB)] 
+kB[{nB + l)Piioi(nA,nB + 1) - nsPnoi (nA, n^)] 
+gti[Piim{nA -I^ub) - PiiminA^riB)] 
+9oi [-Piioi {nA, nB - I) - Pnoi (?^A, (S9) 



dPiooi{nA,nB) _ 
dt ~ 

--^[nA{nA - l)]Piooi{nA,nB) + /aa-Pqooi - 2,nB) 
+-^[("'5 + 2)(ns + l)]Puoi{nA,nB + 2) - /ab-Plooi("-a, "-b) 
+ 2)(nA + l)]Pioii(nA + 2,71^) - /ba-Plooi("-a, "-s) 

^B B 

^[nB{nB - l)]Piooi{nA,nB) + /bb-Piooo(?M, "-s - 2) 

+^'a[(?^A + 1)^1001 (f^A + - ?^A-PlOOl(^^A,^^B)] 

+A:B[(nB + l)Plool(^^A,^^B + 1) - "■B-pL00l(^^A, 

+gio[Piooi{nA -l,nB) - -Piooi("'A, ns)] 
+5M[-Piooi(f^A,f^B - 1) - Piooi{nA,nB)] (SIO) 
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dPoioi{nA,nB] 



+-^[(nyi + 2){nA + l)]Piioi (nA + 2,?!^) - /aa-Poioi("'A, ^^b) 

--^[iT'B{nB - l)]Poioiin A, ub) + fABPoooi{nA,nB - 2) 

+ 2)(nA + l)]PouiinA + 2,725) - /ba-Poioi("'A, 

- l)]Pmm{nA,nB) + fBBPoioo{nA,nB - 2) 

+kA[{nA + l)Poioi(f^A + l,'"'^) - nAPoioiinA,nB)] 
+kB[{nB + l)PolOl(^^A,^^B + 1) - ?^B-PolOl(^^A, 

+5^^! [Poioi (f^A - - Poioi("'A,?^b)] 

+goi[Pmoi{nA,nB - 1) - -PolOl(^^A, (Sll) 

(iPoooi(nA,ns) 



dt 

+-y^[("-A + 2)(nA + l)]Piooi("-A + 2,725) - /aa-Poooi("-a,?t-b) 
+-^[{i^B + 2)(7i5 + l)]PolOl("'A,?^B + 2) - /ab-Poooi("-a,"-b) 
+-|^[("A + 2)(nA + l)]Pooii('^A + 2,715) - fBAPmQi{nA,nB) 

^B B 

^[nB{nB - l)]-Poooi('^A,ra5) + fBBPoooo{nA,nB - 2) 

+kA[{nA + 1)^0001 (f^A + - ?^A-foool(^^A,^^B)] 

+A:B[(n5 + l)Poooi(f^A,f^B + 1) - "■B-Poooi(?^A, ^^B)] 
+500 [^0001 (f^A - l,n5) - -Poooi("'A,n5)] 
+9oi[Poooi{nA,nB - 1) - -Poooi(?^A, ?^b)] (S12) 
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dPiioo{nA,nB) 



dt 

hAA 



2 

hAB 



[uAiriA - l)]Piioo("-A, ns) + /AA-PolOo(^^A - 2, ub) 
nsinB - l)]Plloo(^^A,?^B) + fABPioooinA,nB - 2) 



+ + 2)(n^ + l)]Piiio(nA + 2,nB) - /ba-Piioo(?M, 

-[{uB + 2){nB + l)]PiiQi{nA,nB + 2) - /bbPiioo("-a, 



2 

+fe^[(nA + l)Piioo(nA + l,"-^) - nAPiim{nA,nB)] 
+kB[{nB + l)Piioo('^A,'^B + 1) - nBPiwo{nA,nB)] 
+gu[Pnoo{nA - l,nB) - Piwo{nA,nB)] 
+9oo[Piioo{nA,nB - 1) - Piwo{nA,nB)] (S13) 



dPiooo{nA,nB) _ 
dt ~ 

--^[^A{nA - 1)]Piooo("-a,?t-b) + fAAPoQOoinA - 2,nB) 
+ 2)(nB + l)]Puoo{nA,nB + 2) - /ab-Piooo(?M, 
+ + 2)(n^ + l)]Pioio(raA + 2,nB) - /ba-Piooo(?M, "-s) 

^BB 

-[{uB + 2){nB + l)]Piooi{nA,nB + 2) - /bbPiooo(?M, 



2 

+^a[("'A + l)-Piooo(f^A + ^,nB) - nAPioooi^AjnB)] 
+kB[inB + l)Piooo{nA,nB + 1) - nBPiooo(?^A, -raB)] 
+5io[-Piooo("'A - - -Piooo("'A,'raB)] 

+5'TO[-Pl000("-A,raB - 1) - -Pl000(?^A,raB)] (S14) 
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dPoioo{nA,nB) _ 
dt ~ 

+ + 2){nA + l)]Piioo("-A + 2,71^) - /Ayi-Poioo("-yl, "-b) 

--^[n-BinB - l)]PolOo(^^A,?^B) + /AB-Poooo(?^A, - 2) 
+ + 2)(n^ + l)]Poiio('^A + 2,725) - fBAPoioo{nA,7iB) 

-[{uB + 2){nB + l)]PolOl(?^A,?^B + 2) - /bb-Poioo("-a, 



2 

+fe^[(?T.A + l)PoioQ{nA + l,nB) - nAPQim{nA,nB)] 
+kB[{nB + l)Poioo{nA,nB + 1) - nBPoioo{nA,nB)] 
+5oi[-Poioo("'A - l,™^) - Poioo{nA,nB)] 
+9oo[Poioo{nA,nB - 1) - -PolOo(^^A, -raB)] (S15) 



dPoooo{nA,nB) 



dt 

+-^[{iT'A + 2)(n^ + l)]Piooo("-A + 2,715) - fAAPoooo{nA,nB) 
+ 2)(?iB + l)]Poiooin A, riB + 2) - /ab-Poooo(?M, ?^b) 
+ 2)(n^ + l)]Pooio('^A + 2,72b) - fBAPoooo{nA,nB) 

^BB 

-[{uB + 2){nB + l)]Poooi{nA,nB + 2) - /bbPoooo(?M, ?^b) 



2 

+^a[("'A + l)-Poooo(f^A + ^,nB) - nAPoooo{nA,nB)] 
+kB[inB + l)Poooo{nA,nB + 1) - 71bPoooo("'A, riB)] 
+5'oo[-Poooo("'A - 1,™^) - -Poooo(?^A,?^B)] 
+5'a)[-Poooo(nA, ti^ - 1) - Poooo{nA, ^b)] (S16) 
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Figure 1: Network diagram of canonical gene regulatory circuit of two mu- 
tually opposing proteins that positively self-regulate themselves. Two types 
of genes, A and B are translated into proteins A and B respectively. The 
proteins A{B) can bind to the promoter of the gene A{B) to activate the 
synthesis rate of A{B), which makes a self-activation feedback loop. The 
proteins A{B) can bind to the gene B(A) to repress the synthesis rate of 
B(A), which makes a mutual repression loop. Both protein A and protein B 
bind on promoters as a dimer with the binding rate haA = ^haA^A^nA — 1), 
haB = \haBnB{nB — 1) respectively, and the unbinding rate faA, faS re- 
spectively, with a = {A^B). 
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Figure 2: The potential landscape (contour view) in n^-UB plane for differ- 
ent self activation strength Fa and binding/unbinding speed w. Differentia- 
tion happens with the decrease of the binding/unbinding speed oj (from left 
to right) or the decrease of the activation strength (from top to bottom). 
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(a) Fa = 20, o; = 1000 (b) Fa = 20, o; = 1 (c) Fa = 20, = 0.001 




(d) Fa = 13, oj = 1000 (e) Fa = 13, o; = 1 (f ) Fa = 13, o; = 0.001 




(g) Fa = 0, a; = 1000 (h) Fa = 0, = 1 (i) Fa = 0, cj = 0.001 



Figure 3: The potential landscape (3 dimensional view) in nj^-ns plane 
for different self activation strength and binding/unbinding speed co. 
Differentiation happens with the decrease of the binding/unbinding speed uj 
(from left to right) or the decrease of the activation strength Fa (from top 
to bottom). 
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Figure 4: The MFPT of the differentiation and reprogramming for different 
self activation strength Fa and binding/unbinding speed u. 
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Figure 5: Transition paths for differentiation (blue) and reprogramming 
(purple) with k = 0.1, with self activation strength Fa = 20. 
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